Introduction à la modélisation statistique bayésienne

Ladislas Nalborczyk & Thierry Phénix

Planning


Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 & 4: Modèle de régression linéaire

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Modèles multi-niveaux

Cours 9: Introduction à brms

Cours 10: Data Hackaton


Plan

  • Rappels et objectifs
  • Définition du modèle
  • La distribution postérieure
  • Utilisation du modèle

Rappels

PRINCIPES DE L'ANALYSE DE DONNÉES BAYESIENNE

  • Un ensemble de données / mesures observées à analyser :
    \( d_1, d_2, \dots, d_N \)

  • Un modèle génératif défini par un ensemble de paramètres
    \( \theta_1, \theta_2, \dots, \theta_K \)

  • Une connaissance a priori de ces paramètres

Nous avons les observations

The Bayesian approach places uncertainty not in the observations but rather in one’s lack of knowledge. For a Bayesian, the observed data are not uncertain—you observed what you observed.

Feinberg and Gonzalez (2012)

Rappels

PRINCIPES DE L'ANALYSE DE DONNÉES BAYESIENNE

  • Un ensemble de données / mesures observées à analyser :
    \( d_1, d_2, \dots, d_N \)

  • Un modèle génératif défini par un ensemble de paramètres
    \( \theta_1, \theta_2, \dots, \theta_K \)

  • Une connaissance a priori de ces paramètres

Nous cherchons

Rappels

UTILISATION DE LA RÈGLE DE BAYES

Ces observations et ces paramètres sont regroupés dans un même modèle bayésien :

la distribution conjointe \( P(d_{1:N}~\theta_{1:K}) \).

2 possibilités pour décomposer cette distribution à l'aide du théorème de Bayes :

  • En faisant apparaître la fonction de vraissemblance (likelihood) et la distribution a priori (prior)

\[ P(d_{1:N}~\theta_{1:K})~=~\color{orangered}{P(d_{1:N}~|~\theta_{1:K})}~\color{steelblue}{P(\theta_{1:K})} \]


  • En faisant apparaître la distribution postérieure

\[ P(d_{1:N}~\theta_{1:K})~=~\color{purple}{P(\theta_{1:K}~|~d_{1:N})}~\color{green}{P(d_{1:N})} \]


Ces deux décompositions sont égales !!!

Rappels

INFÉRENCE BAYESIENNE

On infére les nouvelles valeurs des paramètres à partir des données

\[ \color{purple}{P(\theta_{1:K}~|~d_{1:N})} = \frac{\color{orangered}{P(d_{1:N}~|~\theta_{1:K})}~\color{steelblue}{P(\theta_{1:K})}}{\color{green}{P(d_{1:N})}} \]

Ou encore

\[ \color{purple}{P(\theta_{1:K}~|~d_{1:N})} \propto \color{orangered}{P(d_{1:N}~|~\theta_{1:K})}~\color{steelblue}{P(\theta_{1:K})} \]

L'objectif de l'analyse de données bayésienne

Passer d'une connaissance a priori sur les paramètres \( \color{steelblue}{P(\theta_{1:K})} \)
À une connaissance à postériori \( \color{purple}{P(\theta_{1:K}~|~d_{1:N})} \)

Rappels

ANALYSE DE DONNÉES BAYÉSIENNE

1. Définir le modèle

  • Identifier les paramètres du modèle génératif (Likelihood) associés aux données
  • Définir une distribution a priori (prior) sur ces paramètres


2. Mettre à jour le modèle

  • Calculer la distribution a posteriori des paramètres (ou une bonne approximation)


3. Interpréter la distribution postérieure (répondre à la question posée)

  • Estimation du paramètre
  • Comparaison de modèles


4. Vérification prédictive du modèle

  • S'assurer que le modèle génère des données cohérentes avec les données observées

Objectifs de ce cours





Illustrer les différents aspects de cette démarche

à l'aide d'un modèle simple


le modèle Beta-Binomial


Planning


Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 & 4: Modèle de régression linéaire

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Modèles multi-niveaux

Cours 9: Introduction à brms

Cours 10: Data Hackaton


Plan

  • Rappels et objectifs
  • Définition du modèle
    • Likelihood
    • Prior
  • La distribution postérieure
  • Utilisation du modèle

Le modèle beta-binomial

POURQUOI COMMENCER PAR CE MODÈLE

  • Le modèle beta-binomial couvre un grand nombre de problèmes de la vie courante :

    • Réussite / Echec à un test psychométrique ?
    • Présence / absence d'effets secondaires lors d'un test d'un médicament ?
    • Résultat d'un questionnaire à réponse binaire vrai/faux
    • Estimation des résultats du deuxième tour de l'éléction présidentielle
  • C'est un modèle simple

    • Un unique paramètre
    • Solution analytique

Le modèle beta-binomial

FIL ROUGE - lancer d'une pièce de monnaie

Définition du modèle - Les données

  • Génération des données, la valeur de \( \theta \) étant connue (ici, \( \theta=0.6 \)).
library(tidyverse)
set.seed(1234567)
sample(x = c(0, 1), size = 50, replace = TRUE, prob = c(0.4,0.6)) %>%
        data.frame %>%
        mutate(x = seq_along(.), y = cumsum(.) / seq_along(.) ) %>%
        ggplot(aes(x = x, y = y), log = "y") +  geom_line(lwd = 1, col = "steelBlue") +
        geom_hline(yintercept = 0.5, lty = 3) +
        xlab("nombre de lancers") + ylab("proportion de faces") +
        ylim(0, 1) + theme_bw(base_size = 20)

plot of chunk unnamed-chunk-1

Définition du modèle

LIKELIHOOD

Modèle d'un lancer unique

Jacques Bernoulli

Jacques ou Jakob Bernoulli (1654-1705)
mathématicien et physicien suisse
(né et mort à Bâle)

Modèle d'un lancer unique

Depuis Bernoulli, on sait calculer la probabilité du résultat d'un lancer, du moment que l'on connait le biais de la pièce

\[ P(Y=y~|~\theta)=\theta^y~(1-\theta)^{(1-y)} \]

avec \( \theta \) la probabilité de l'évènement “le résultat du lancer est \( \textit{Face} \)”


Modèle d'une suite de lancers (processus de Bernoulli)

Pour une série de lancers \( \{Y_i\} \) indépendants et identiquement distribués (i.i.d.)

\[ P(\{Y_i=y_i\}~|~\theta)=\theta^z~(1-\theta)^{(N-z)} \]

avec \( z=\sum_i y_i \) (nombre de succès) et N nombre d'observations

Définition du modèle

LIKELIHOOD

  • Nous considérons \( z \) comme étant le nombre de succès \( \{y_i\} \)
  • Nous considérons N (le nombre d'observations) comme étant une constante
  • Nous considérons \( \theta \) comme étant un paramètre de notre modèle


La fonction de likelihood s'écrit :

\[ L(\theta~|~z~N)=P(z~|~\theta~N)=\theta^z~(1-\theta)^{(N-z)} \]

c'est une fonction de la variable \( \theta \)

Définition du modèle

LIKELIHOOD - distribution Binomiale

Exemple :
Représentation graphique de la fonction de vraissemblance de \( \theta \) pour 10 lancers et 7 Faces

NbrTrial     = 100
N            = 10
nbrSuccess   = 7
theta        = seq(from = 0, to = 1, by = 1/NbrTrial)
likelihood   = (theta^(nbrSuccess)) * (1 - theta)^(N - nbrSuccess) 
data.frame(theta, likelihood) %>%
ggplot(aes(x = theta, y = likelihood)) + 
    geom_area(color="orangered", fill = "orangered", alpha =.4) +
    theme_bw(base_size = 20) + theme(legend.position="none",
          text = element_text(size=30))

plot of chunk unnamed-chunk-2

Définition du modèle

PRIOR

Comment définir un prior dans le cas du lancer de pièce ?

Aspect sémantique, rendre compte ou modéliser :

  • Une absence d'information
  • La connaissance d'observations antérieures concernant la pièce lancée.
  • Un niveau de certitude concernant ces observations antérieures


Aspect mathématique, pour une solution entièrement analytique :

  • Les distribution postérieure et prior doivent avoir la même forme
  • La vraissemblance marginale doit se calculer analytiquement

Définition du modèle

PRIOR - distribution Beta

Un candidat : La distribution beta

La distribution beta:

\[ \begin{align} p(\theta~|~a,b) &= beta(\theta~|~a,b) \\ &= \theta^{(a-1)}(1-\theta)^{(b-1)}/B(a,b) \\ &\propto \theta^{(a-1)}(1-\theta)^{(b-1)} \end{align} \]

où \( a \) et \( b \) sont deux paramètres tels que \( a\geq0~et~b\geq0 \) et \( B(a,b) \) est une constante de normalisation


Attention : Ne pas confondre dans R : la fonction de \( \theta \) et la constante (indépendante de \( \theta \))

a     = 2
b     = 3
theta = seq( from= 0 , to= 1 , length.out= 7 )
dbeta(theta,a,b)
[1] 0.0000000 1.3888889 1.7777778 1.5000000 0.8888889 0.2777778 0.0000000
beta(a,b)
[1] 0.08333333

Définition du modèle

PRIOR - distribution Beta

Vérifions les contraintes mathématiques

Le résultat du produit d'une distribution beta et d'une fonction de vraissemblance de Bernoulli est proportionnel à une distribution beta.
On dit que la distribution beta est la distribution conjugée de la fonction de vraissemblance de Bernoulli.


Soit un prior définit par \( ~~~\color{steelblue}{p(\theta~|~a,b) \propto \theta^{(a-1)}(1-\theta)^{(b-1)}} \)
Soit une fonction de vraissemblance associée à z “Face” pour N lancers \( ~~\color{orangered}{P(\{y_i\}~|~\theta)=\theta^z~(1-\theta)^{(N-z)}} \)

\[ \begin{align} \color{purple}{p(\theta~|~\{y_i\})} &\propto \color{orangered}{P(\{y_i\}~|~\theta)}~\color{steelblue}{p(\theta~|~a,b)} &&\mbox{inférence Bayésienne} \\ &\propto \color{orangered}{\theta^z~(1-\theta)^{(N-z)}}~\color{steelblue}{\theta^{(a-1)}(1-\theta)^{(b-1)}} &&\mbox{application des formules précédentes} \\ &\propto \color{purple}{\theta^{(z+a-1)}~(1-\theta)^{(N-z+b-1)}} &&\mbox{en regroupant les termes identiques} \\ &\propto \color{purple}{\theta^{(A-1)}(1-\theta)^{(B-1)}} &&\mbox{avec } A=z+a \mbox{ et } B=N-z+b \end{align} \]

Conclusions : la distribution a priori et la distribution a posteriori ont la même forme

Définition du modèle

PRIOR - distribution Beta

Vérifions les aspects sémantiques

  • On peut exprimer l'absence de connaissance à priori - a=1 et b=1 (courbe bleue).
  • On modélise un biais en faveur de Face - a > b (courbe orange).
  • On modélise un biais en faveur de Pile - a < b (courbe verte).


plot of chunk unnamed-chunk-4

Définition du modèle

PRIOR - distribution Beta

Vérifions les aspects sémantiques

  • Le niveau de certitude augmente avec la somme (\( \kappa = a+b \))

    • Aucune idée sur la provenance de la pièce : a = b = 1 -> prior faible (courbe bleue)
    • En attendant le début de l'expérience, on a lancé la pièce 10 fois pour 5 “Face” : a = b = 5 (courbe orange)
    • La pièce provient de la banque de France : a = b = 50 -> prior fort (courbe verte)

plot of chunk unnamed-chunk-5

Définition du modèle

PRIOR - distribution Beta

Le plus souvent, on dispose :

  • D'information sur la tendance centrale du paramètre (mode, moyenne)
  • D'un niveau de certitude ( \( \kappa = a+b \) )

Définition du modèle

PRIOR - distribution Beta

Le mode
Supposons que l'on ait une estimation de la valeur la plus probable \( \omega \) du paramètre \( \theta \)
On peut définir les paramètres a et b de notre prior

\[ \begin{align} a &= \omega(\kappa-2)+1 \\ b &= (1-\omega)(\kappa-2)+1 &&\mbox{pour } \kappa>2 \end{align} \]


Pour \( ~~\omega=0.65 \) et \( ~~\kappa=25 \) alors \( ~~~p(\theta)=beta(\theta~|~15.95, 9.05)~ \) (courbe orange)
Pour \( ~~\omega=0.65 \) et \( ~~\kappa=10 \) alors \( ~~~p(\theta)=beta(\theta~|~6.2, 3.8)~ \) (courbe bleue)

plot of chunk unnamed-chunk-6

Définition du modèle

PRIOR - distribution Beta

La moyenne
Supposons que l'on ait une estimation de la valeur moyenne \( \mu \) du paramètre \( \theta \)
On peut définir les paramètres a et b de notre prior

\[ \begin{align} a &= \mu \kappa \\ b &= (1-\mu) \kappa &&\mbox{ pour } \kappa>2 \end{align} \]


Pour \( \mu=0.8 \) et \( \kappa=7 \) alors \( p(\theta)=beta(\theta~|~5.6,1.4) \) (courbe orange)
Pour \( \mu=0.8 \) et \( \kappa=22 \) alors \( p(\theta)=beta(\theta~|~17.6,4.4) \) (courbe bleue)

plot of chunk unnamed-chunk-7

Définition du modèle

PRIOR - distribution Beta

La moyenne et la variance
Supposons que l'on ait une estimation a priori de la valeur moyenne \( \mu \) et de la variance \( \sigma \) du paramètre \( \theta \)
On peut définir les paramètres a et b de notre prior

\[ \begin{align} a &= \mu~(\mu~(1-\mu)~/~\sigma^2-1) \\ b &= (1-\mu)~\mu~(\mu~(1-\mu)~/~\sigma^2-1) ,~~~~~ \mbox{ pour } \sigma < 0.28867 \\ \end{align} \]


Pour \( ~~\mu=0.8 \) et \( ~~\sigma=0.15 \) alors \( ~~~p(\theta)=beta(\theta~|~5.92, 3.19) \) (orange)
Pour \( ~~\mu=0.8 \) et \( ~~\sigma=0.1 \) alors \( ~~~p(\theta)=beta(\theta~|~14.14, 7.61) \) (bleue)

plot of chunk unnamed-chunk-8

Définition du modèle

PRIOR - Son influence sur la moyenne de la distribution postérieure

Pour une distribution a priori \( \color{steelblue}{beta(\theta~|~a,b)} \), de moyenne \( \color{steelblue}{\frac{a}{a+b}} \), et un résultat de \( \color{orangered}{z} \) Face sur \( \color{orangered}{N} \) lancers, la moyenne de la distribution postérieure est :

\[ \frac{z+a}{N+a+b} = \color{orangered}{\frac{z}{N}}~\frac{N}{N+a+b}+\color{steelblue}{\frac{a}{a+b}}~\frac{a+b}{N+a+b} \]

Cas \( \mbox{ } N < a+b \)

plot of chunk unnamed-chunk-9

Définition du modèle

PRIOR - Son influence sur la moyenne de la distribution postérieure

Pour une distribution a priori \( \color{steelblue}{beta(\theta~|~a,b)} \), de moyenne \( \color{steelblue}{\frac{a}{a+b}} \), et un résultat de \( \color{orangered}{z} \) Face sur \( \color{orangered}{N} \) lancers, la moyenne de la distribution postérieure est :

\[ \frac{z+a}{N+a+b} = \color{orangered}{\frac{z}{N}}~\frac{N}{N+a+b}+\color{steelblue}{\frac{a}{a+b}}~\frac{a+b}{N+a+b} \]

Cas \( \mbox{ } N = a+b \)

plot of chunk unnamed-chunk-10

Définition du modèle

PRIOR - Son influence sur la moyenne de la distribution postérieure

Pour une distribution a priori \( \color{steelblue}{beta(\theta~|~a,b)} \), de moyenne \( \color{steelblue}{\frac{a}{a+b}} \), et un résultat de \( \color{orangered}{z} \) Face sur \( \color{orangered}{N} \) lancers, la moyenne de la distribution postérieure est :

\[ \frac{z+a}{N+a+b} = \color{orangered}{\frac{z}{N}}~\frac{N}{N+a+b}+\color{steelblue}{\frac{a}{a+b}}~\frac{a+b}{N+a+b} \]

Cas \( \mbox{ } N > a+b \)

plot of chunk unnamed-chunk-11

Définition du modèle

CONCLUSION

There are two main desiderata for a mathematical description of data.

  • First, the mathematical form should be comprehensible with meaningful parameters.
  • The second desideratum for a mathematical description is that it should be descriptively adequate, which means, loosely, that the mathematical form should “look like” the data.

Kruschke (2015)

Planning


Cours

Cours 1 : Introduction à l'inférence bayésienne

Cours 2 : Modèle beta-binomial

Cours 3 & 4 : Modèle de régression linéaire

Cours 5 : Markov Chain Monte Carlo

Cours 6 : Modèle linéaire généralisé

Cours 7 : Comparaison de modèles

Cours 8 : Modèles multi-niveaux

Cours 9 : Introduction à brms

Cours 10 : Data Hackaton


Plan

  • Rappels et objectifs
  • Définition du modèle
  • La distribution postérieure
    • Cas analytique
    • Méthode Grid
    • Le modèle binomial
  • Utilisation du modèle

La distribution postérieure

INTRODUCTION

The posterior distribution is always a compromise between the prior distribution and the likelihood function.
Kruschke (2015)

plot of chunk unnamed-chunk-12

La distribution postérieure

PRINCIPE

\[ Posterior = \frac{Likelihood~\times~Prior}{\color{orangered}{Marginal~Likelhood}} \]


PROBLÈME :

le calcul de la “marginal likelihood”, \( \color{green}{~p(D)} \)

\[ \begin{align} \color{green}{p(D)} &= \int P(D,\theta) \, \mathrm d\theta &&\mbox{marginalisation sur le paramètre} \theta \\ \color{green}{p(D)} &= \color{green}{\int P(D~|~\theta)~P(\theta) \, \mathrm d\theta} && \mbox{Application de la règle du produit} \end{align} \]

Dans la littérature, on le trouve sous le nom de : “evidence”, “probabilité des données” ou “Marginal Likelihood”

La distribution postérieure

PRINCIPE

\[ Posterior = \frac{Likelihood~\times~Prior}{\color{orangered}{Marginal~Likelhood}} \]


Trois méthodes pour résoudre ce problème :


  1. Solution analytique \( ~~\longrightarrow~~ \) Utilisation d'une likelihood conjuguée du prior
  2. Solution discrètisée \( ~~\longrightarrow~~ \) Calcul de la solution sur un ensemble fini de points (Méthode Grid)
  3. Solution approchée \( ~~\longrightarrow~~ \) On échantillonne “intelligemment” sur l'espace conjoint des paramètres (Méthode MCMC)

La distribution postérieure

CAS ANALYTIQUE

Utilisation d'une likelihood conjuguée du prior : cas Beta-binomiale


\[ \begin{align} \color{purple}{p(\theta~|~z,N)} &= \color{orangered}{P(z~|~\theta,N)}~\color{steelblue}{p(\theta)/\color{green}{P(z~|~N)}} &&\mbox{Équation de Bayes} \\ &= \color{orangered}{\theta^z~(1-\theta)^{(N-z)}}~\color{steelblue}{\frac{\theta^{(a-1)}(1-\theta)^{(b-1)}}{B(a,b)}}/\color{green}{P(z~|~N)} &&\mbox{Application des formules précédetes} \\ &= \frac{\color{orangered}{\theta^z~(1-\theta)^{(N-z)}}~\color{steelblue}{\theta^{(a-1)}(1-\theta)^{(b-1)}}}{\color{steelblue}{B(a,b)}\color{green}{P(z~|~N)}} \\ &= \frac{\theta^{(\color{orangered}{z}+\color{steelblue}{a-1})}~(1-\theta)^{(\color{orangered}{N-z}+\color{steelblue}{b-1})}}{\color{steelblue}{B(a,b)}\color{green}{P(z~|~N)}} &&\mbox{En regroupant les termes identiques} \\ &= \frac{\theta^{(\color{orangered}{z}+\color{steelblue}{a-1})}~(1-\theta)^{(\color{orangered}{N-z}+\color{steelblue}{b-1})}}{B(\color{orangered}{z}+\color{steelblue}{a},\color{orangered}{N-z}+\color{steelblue}{b})} &&\mbox{Par identification à une fonction beta} \\ \end{align} \]


L'avantage d'avoir un prior conjugué nous dispense du calcul du dénominateur

La distribution postérieure

CAS ANALYTIQUE

Distributions discrètes

La distribution postérieure

CAS ANALYTIQUE

Distributions continues


TROP CONTRAIGNANT !!!

Le modèle (likelihood + prior) doit être défini à partir de l'interprétation que l'on peut faire des paramètres et des distributions
et non pour faciliter les calculs

La distribution postérieure

MÉTHODE GRID

  1. Définir la grille
  2. Calculer la valeur du prior pour chaque valeur de la grille
  3. Calculer la valeur de la likelihood pour chaque valeur de la grille
  4. Calculer le produit prior x likelihood pour chaque valeur de la grille, puis normalisation du résultat


plot of chunk unnamed-chunk-13

La distribution postérieure

MÉTHODE GRID

  1. Définir la grille
  2. Calculer la valeur du prior pour chaque valeur de la grille
  3. Calculer la valeur de la likelihood pour chaque valeur de la grille
  4. Calculer le produit prior x likelihood pour chaque valeur de la grille, puis normalisation du résultat


plot of chunk unnamed-chunk-14

La distribution postérieure

MÉTHODE GRID

  1. Définir la grille
  2. Calculer la valeur du prior pour chaque valeur de la grille
  3. Calculer la valeur de la likelihood pour chaque valeur de la grille
  4. Calculer le produit prior x likelihood pour chaque valeur de la grille, puis normalisation du résultat


plot of chunk unnamed-chunk-15

La distribution postérieure

MÉTHODE GRID

  1. Définir la grille
  2. Calculer la valeur du prior pour chaque valeur de la grille
  3. Calculer la valeur de la likelihood pour chaque valeur de la grille
  4. Calculer le produit prior x likelihood pour chaque valeur de la grille, puis normalisation du résultat


plot of chunk unnamed-chunk-16

La distribution postérieure

MÉTHODE GRID - grille affinée

  1. Définir la grille
  2. Calculer la valeur du prior pour chaque valeur de la grille
  3. Calculer la valeur de la likelihood pour chaque valeur de la grille
  4. Calculer le produit prior x likelihood pour chaque valeur de la grille, puis normalisation du résultat


plot of chunk unnamed-chunk-17

La distribution postérieure

MÉTHODE GRID - Problème du nombre de paramètres


En affinant la grille on augmente le temps de calcul

  • 3 paramètres avec une grille de \( 10^3 \) noeuds => une grille de \( 10^9 \) points de calcul
  • 10 paramètres avec une grille de \( 10^3 \) noeuds => une grille de \( 10^{30} \) points de calcul

Le “superordinateur” chinois Tianhe-2 réalise \( 33,8 x 10^{15} \) opérations par seconde
Si on considère qu'il réalise 1 opération par noeud de la grille, il lui faudrait \( 2.96 x 10^{13} \) pour parcourir la grille une fois !!!
Par comparaison, l'âge de l'univers est approximativement de \( (4,354 ± 0,012)×10^{17} \) secondes

La distribution postérieure

MÉTHODE PAR ÉCHANTILLONNAGE - Le principe

  • On ne calcule pas explicitement la distribution postérieure
  • On en donne une estimation à partir du produit likelihood x prior calculée sur un échantillon des noeuds de la grille

En pratique :

  • La fonction sample de R
trajLength     = 250 # nbr of measure                    
theta          = 1:7
ptheta         = theta
trajectory     = sample(theta, prob = ptheta, size = trajLength, replace = TRUE)
  • Les méthodes MCMC, il existe différentes implémentations.
    (Metropolis, Gibbs, Hamilton, …)

La distribution postérieure

MÉTHODE PAR ÉCHANTILLONNAGE - precision

plot of chunk unnamed-chunk-19

plot of chunk unnamed-chunk-20

La précision dépend de la taille de l'échantillon

La distribution postérieure

EN PRATIQUE

  • Cas analytique : Une fonction
p_grid    =  seq( from=0 , to=1 , length.out=1000 )
a = b     = 1
N         = 9
z         = 6
posterior = dbeta(p_grid, z+a-1, N-z+b-1)
  • Cas Discret : Une valeur par noeud de la grille
p_grid      = seq( from=0 , to=1 , length.out = 1000 )
prior       = rep( 1 , 1000 )
likelihood  = dbinom( 6 , size=8 , prob=p_grid )
posterior   = likelihood * prior
posterior   = posterior / sum(posterior)
  • Cas échantillons : Un ensemble de valeurs par noeud de la grille
samples = sample( p_grid , prob = posterior , size = 1e4 , replace = TRUE )

La distribution postérieure

CONCLUSION

Méthode analytique

  • La distribution postérieure est décrite explicitement
  • Le modèle est fortement contraint

Méthode Grid

  • La distribution postérieure n'est donnée que pour un ensemble finie de valeurs
  • Plus la grille est fine, meilleure est l'estimation de la distribution postérieure
  • Compromis Précision - Temps de calcul => Faible nombre de paramètres

Méthodes par échantillonnage (MCMC)

  • La distribution postérieure est décrite par un échantillon
  • Vérification résultat

La distribution postérieure est presque toujours représentée par un échantillon

Planning


Cours

Cours 1 : Introduction à l'inférence bayésienne

Cours 2 : Modèle beta-binomial

Cours 3 & 4 : Modèle de régression linéaire

Cours 5 : Markov Chain Monte Carlo

Cours 6 : Modèle linéaire généralisé

Cours 7 : Comparaison de modèles

Cours 8 : Modèles multi-niveaux

Cours 9 : Introduction à brms

Cours 10 : Data Hackaton


Plan

  • Rappels et objectifs
  • Définition du modèle
  • La distribution postérieure
  • Utilisation du modèle
    • Estimation des tendances centrales
    • Highest density interval (HDI)
    • ROPE
    • Bayes factor
    • BIC

Utilisation du modèle

Bayesian analysis Tendancies

Utilisation du modèle

ESTIMATION DE LA TENDANCE CENTRALE

DESCRIPTION : A partir d'un échantillon d'une distribution postérieure, on peut calculer la moyenne, le mode et la médiane

library(rethinking)
chainmode( samples , adj=0.01 ) # pour le mode ou MAP
mean( samples )                 # pour la moyenne
median( samples )               # pour la médiane

Exemple pour un prior uniforme, et 10 lancers successifs et 3 Faces plot of chunk unnamed-chunk-25

Utilisation du modèle

ESTIMATION DE LA TENDANCE CENTRALE - Cas “pathologique”“

DESCRIPTION : A partir d'un échantillon d'une distribution postérieure, on peut calculer la moyenne, le mode et la médiane

library(rethinking)
chainmode( samples , adj=0.01 ) # pour le mode ou MAP
mean( samples )                 # pour la moyenne
median( samples )               # pour la médiane

Exemple pour un prior uniforme, et 3 lancers successifs et 3 Faces plot of chunk unnamed-chunk-27

Utilisation du modèle

ESTIMATION DE LA TENDANCE CENTRALE

DÉCISION : Comment chosir la valeur qui résume le “mieux” la distribution postérieure
Utiliser des fonctions d'espérance de gain ou, ce qui revient au même, de perte attendue

# pour la moyenne
sapply( p_grid , function(d) sum( posterior*( d - p_grid )^2 ) )
# pour la médiane
sapply( p_grid , function(d) sum( posterior*abs( d - p_grid ) ) )


plot of chunk unnamed-chunk-29

Utilisation du modèle

HIGHEST POSTERIOR DENSITY INTERVAL (HPDI)

  • Le HPDI indique les valeurs du paramètre qui sont le plus crédibles
  • Plus le HPDI est étroit et plus le degré de certitude est élevé
  • La largeur du HPDI diminue avec l'augmentation du nombre de mesures
  • Le pourcentage de la masse de probabilité n'est pas nécessairement 95


Définition: les valeurs du paramètre x contenues dans le HDI 89% sont telles que \( p(x)>W \) et

\[ \int_{x:p(x)>W} p(x) \, \mathrm dx=0.89 \]

Utilisation du modèle

HIGHEST POSTERIOR DENSITY INTERVAL (HPDI)

Bayesian analysis Tendancies

Utilisation du modèle

HIGHEST POSTERIOR DENSITY INTERVAL (HPDI)

En pratique - Package BEST


library(coda)
library(BEST)
set.seed(1234567)
iterMax = 10000
p_grid  = seq( from=0 , to=1 , length.out=iterMax )
pTheta  = dbeta( p_grid, 3, 10)
massVec = pTheta/sum(pTheta)
samples = sample( p_grid , size=1e4 , replace=TRUE , prob=pTheta )

#plotPost(samples, credMass = 0.85)
hdi( samples , credMass = 0.89 )
     lower      upper 
0.05280528 0.39613961 
attr(,"credMass")
[1] 0.89

Utilisation du modèle

HIGHEST POSTERIOR DENSITY INTERVAL (HPDI)

En pratique - Package BEST


plot of chunk unnamed-chunk-31

Utilisation du modèle

HIGHEST POSTERIOR DENSITY INTERVAL (HPDI)

En pratique - Package rethinking

library(rethinking)
set.seed(1234567)
iterMax = 10000
p_grid  = seq( from=0 , to=1 , length.out=iterMax )
pTheta  = dbeta( p_grid, 3, 10)
massVec = pTheta / sum(pTheta)
samples = sample( p_grid , size = 1e4 , replace = TRUE , prob = pTheta )

# utilise le package CODA pour calculer HPDI
HPDI( samples, prob = 0.89)
     |0.89      0.89| 
0.05280528 0.39613961 

Utilisation du modèle

ROPE - Region Of Practical Equivalence

On l'utilise pour “tester une hypothèse”


  • Pour tester H0, la décision est basée sur le HPDI et sur le ROPE
  • La valeur du paramètre (ici 0) est rejetée si son ROPE exclus le 95% HPDI
  • La valeur du paramètre (ici 0) est accepté si le 95% HDI est entièrement inclus dans son ROPE
  • Si le HDI et le ROPE se chevauchent on ne peut rien dire…

Utilisation du modèle

ROPE - Region Of Practical Equivalence


plot of chunk unnamed-chunk-33

Utilisation du modèle

COMPARAISON DE MODÈLES - Principes

On considère une pièce. On la lance 200 fois et on obtient 115 Faces.

QUESTION : La pièce est-elle biaisée?


Nous construisons deux modèles et essayons de savoir lequel est le plus probable

\[ \begin{eqnarray*} \left\{ \begin{array}{ll} M_0~:~X \sim \mathcal{B}(N,0.5) & \quad \text{la pièce n'est pas biaisée}\\ M_1~:~X \sim \mathcal{B}(N,\theta) & \quad \text{la pièce est biaisée} \end{array}\right. \end{eqnarray*} \]

Utilisation du modèle

COMPARAISON DE MODÈLES - Bayes Factor

Le facteur de Bayes fait le rapport des vraissemblances des deux modèles.

\[ \frac{P(M_0~|~D)}{P(M_1~|~D)}=\color{orangered}{\frac{P(D~|~M_0)}{P(D~|~M_1)}}~\frac{P(M_0)}{P(M_1)} \]


Dans cet exemple, le calcul de \( P(D~|~M_0) \) est facile

\( P(D~|~M_0) = C_{200}^{115}~(0.5)^{115}~(1-0.5)^{(200-115)} = 0.005955 \)

Utilisation du modèle

COMPARAISON DE MODÈLES - Bayes Factor

Le facteur de Bayes fait le rapport des vraissemblances des deux modèles.

\[ \frac{P(M_0~|~D)}{P(M_1~|~D)}=\frac{P(D~|~M_0)}{\color{orangered}{P(D~|~M_1)}}~\frac{P(M_0)}{P(M_1)} \]


Pour calculer de \( P(D~|~M_1) \), c'est plus compliqué : \( \color{orangered}{~~~\textbf{On ne connait pas }\theta} \)

Le calcul de la vraissemblance dans ce cas se fait par marginalisation sur le paramètre \( \theta \) \[ \begin{align} P(D~|~M_1) &= \int_0^1 P(D~|~M_1, \theta) P(\theta) \, \mathrm d\theta \\ &= \int_0^1 C^z_N~\theta^z~(1-\theta)^{N-z}~\text{d}\theta \\ &= N+1 \end{align} \]

Utilisation du modèle

COMPARAISON DE MODÈLES - Bayes Factor

Le facteur de Bayes fait le rapport des vraissemblances des deux modèles.

\[ \frac{P(M_0~|~D)}{P(M_1~|~D)}=\color{orangered}{\frac{P(D~|~M_0)}{P(D~|~M_1)}}~\frac{P(M_0)}{P(M_1)} \]


Dans l'exemple,

\[ B_{01} = \frac{P(M_0~|~D)}{P(M_1~|~D)} = \frac{0.005955}{0.005} = 1.1971 \]


Le rapport de probabilités a augmenté de 20% après avoir pris connaissance des données.
La balance penche ici en faveur du modèle nul, \( B_{01} > 1 \)

Utilisation du modèle

COMPARAISON DE MODÈLES - Autre exemple

BIC critère d’information bayésien

Schwarz (1977) propose une approximation du log de la vraissemblance:

La vraissemblance du modèle pas toujours calculable : \[ P(D~|~M)=\int_{\theta} P(D~|~M, \theta) \, \mathrm d\theta \]

Calcul d'une valeur approchée \[ ln(P(D~|~\theta)) \approx ln(L(\widehat{\theta})) - \frac{t}{2}~ln(N) \]

avec \( L(\widehat{\theta}) \) le maximum de vraissemblance
\( t \) le nombre de paramètres
\( N \) le nombre d'observations

Utilisation du modèle

MODEL CHECKING

Les deux rôles de la fonction de vraissemblance

  • C'est une fonction de \( \theta \) pour le calcul de la distribution postérieure : \( L(\theta~|~z~N) \)
  • Lorsque \( \theta \) est fixé, c'est une distribution de probablilité : \( P(z~|~\theta~N)=\theta^z~(1-\theta)^{(N-z)} \)

On peut utiliser cette distribution de probabilité pour générer des données

Exemple
Générer 10000 valeurs à partir d'une loi binomiale basée sur 9 lancers et une probabilité de Face de 0,6 :

w <- rbinom(n = 10e4, size = 9, prob = 0.6)

Utilisation du modèle

MODEL CHECKING

Deux sources d'incertitude dans ces prédictions

  • Incertitude lié à la valeur de \( p \) (évidemment !)
    -> Chaque valeur apparaît avec une probabilité \( p \)
  • Incertitude sur la valeur de \( p \) elle-même
    -> Pour chaque valeur de \( p \) on peut calculer une distribution implicite

Exemple
Générer 10000 valeurs à partir d'une loi binomiale basé sur 9 lancers et une probabilité de Face de 0,6
pour chaque valeur de p

w <- rbinom(n = 10e4, size = 9, prob = samples)

Utilisation du modèle

MODEL CHECKING